#####
#####
###     Run circular regression reported in Table 3
###     and use results to create Figure 5
#####
#####


library(foreign)
library(circular)
library(lattice)

###
###   Read data from Stata dataset,
###   "circreg.dta"
###

circreg <- read.dta(file.choose())

###
###   Calculate angle of each row vector, in radians,
###   using 3:00 position as origin, going counter-
###   clockwise. Also declare the angle to be a 
###   circular object, and create new version of
###   data frame that omits observations with
###   missing values
###

circreg$ang <- atan(circreg$dim2/circreg$dim1) + ((circreg$dim1 < 0) * pi) +
  ((circreg$dim2 <0 & circreg$dim1 > 0) * 2 * pi)

circreg$angle <- circular(circreg$ang)

circreg2 <- na.omit(circreg)

###
###   Recode ideology to seven-point range, from -3 to +3
###

circreg2$ideol <- (circreg2$ideolpst * (6/100)) - 3

###
###   Recode party ID to range from -3 to 3
###

circreg2$postpid7 <- circreg2$postpid7 - 4

###
###   Center the age variable at the mean age in the sample
###

circreg2$age2 <- circreg2$age - mean(circreg2$age)

###
###   Carry out circular regression, producing
###   results reported in Table 3
###

creg1 <- lm.circular(y = circreg2$angle,
                     x = cbind(circreg2$age2, 
                               circreg2$somecollege, circreg2$collegegrad,
                               circreg2$commit,
                               circreg2$postpid7, circreg2$ideol),
                     type = "c-l",
                     init = rep(.008, 6),
                     verbose = T
)

creg1

###   Multiply range of each variable by its coefficient
###   to estimate size of effect
###

(max(circreg2$age) - min(circreg2$age)) * creg1$coefficients[1]

(max(circreg2$commit) - min(circreg2$commit)) * creg1$coefficients[4]

(max(circreg2$postpid7) - min(circreg2$postpid7)) * creg1$coefficients[5]

(max(circreg2$ideol) - min(circreg2$ideol)) * creg1$coefficients[6]

###
###   Calculate vectors for the most extreme
###   segments of the respondents
###

###
###   Create "Most conservative" vector,
###   using results from creg1
###

right <- c(36.47839, 0, 0, 
           12, 3, 3)


###
###   Create "Most liberal" vector
###

left <- c(-30.52161, 1, 1,
          0, -3, -3)

lin.pred.left <- sum(creg1$coefficients * left)

lin.pred.right <- sum(creg1$coefficients * right)

angle.left <- -0.9507 + (2 * atan(lin.pred.left))

angle.left

angle.right <- -0.9507 + (2 * atan(lin.pred.right))

angle.right

###
###   Read value coordinates from file, 
###   "value coords, with names.txt"
###

val.coords <- read.table(file.choose(), header = T)

###
###   obtain bisector of extreme vectors
###

bisectd1 <- (cos(angle.left) + cos(angle.right))/2

bisectd2 <- (sin(angle.left) + sin(angle.right))/2

circreg$side <- "leftward"

circreg$side[circreg$ang > (pi + atan(bisectd2/bisectd1)) & 
               circreg$ang <=  ((2 * pi) + atan(bisectd2/bisectd1))] <- "rightward"

###
###   Calculate mean vectors for subsets of 
###   observations that fall on the left and
###   right sides of the bisector, respectively.
###

meanleftd1 <- mean(circreg$dim1[circreg$side == "leftward"], na.rm = T)

meanleftd2 <- mean(circreg$dim2[circreg$side == "leftward"], na.rm = T)

meanrightd1 <- mean(circreg$dim1[circreg$side == "rightward"], na.rm = T)

meanrightd2 <- mean(circreg$dim2[circreg$side == "rightward"], na.rm = T)

###
###   Plot most extreme vectors, the bisector between the two
###   extremes, and the mean vectors for the subsets of
###   observations that fall on either side of the bisector,
###   along with the value points.
###
###   THIS PRODUCES FIGURE 5
###

xyplot(dim2 ~ dim1, data = circreg,
       aspect = 1,
       panel = function (x, y) {
         panel.xyplot(val.coords$dim1, val.coords$dim2, pch = 16, cex = .75, 
                      col = "black")
         panel.text(val.coords$dim1, val.coords$dim2, 
                    labels = c("Freedom", "Equality", "Economic\nsecurity",  
                               "Morality", "Individualism", "Social\norder",       
                               "Patriotism"),
                    pos = 1, cex = .75)
         panel.segments(rep(0, 2), rep(0, 2), 
                        c(cos(angle.left), cos(angle.right)), c(sin(angle.left), 
                                                                sin(angle.right)), col = "gray")
         panel.abline(a = 0, b = (bisectd2/bisectd1), col = "gray")
         panel.segments(rep(0, 2), rep(0, 2), c(meanleftd1, meanrightd1),
                        c(meanleftd2, meanrightd2), col = "black",
                        lwd = 2)
       },
       xlab = "",
       ylab = "",
       scales = list(draw = FALSE),
       xlim = c(-1, 1),
       ylim = c(-1, 1)
)


